4  CVM ppp

4.1 Pretratamientos de la Información

Cargar Recursos Externos

source("scripts/recursos.R")
source("scripts/funciones.R")

Definición de variables generales

## geographic projections
crs_latlon <- "+proj=longlat +datum=WGS84 +no_defs"
crs_utm <- 
  "+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

Definición de variables específicas

# List of crimes/delitos and cities/ciudades
crimes = c(
  "Minor property offenses",
  "Robberies",
  "Burglaries",
  "Drunkennes, damages and disorders",
  "Family violence and aggression",
  "Injuries, drugs and weapons",
  "High violence and murder"
)
delitos = c("hurto", "robo", "asalto", "amen", "abus", "crim", "viol")

#ciudad de ejemplo
cities <- c("Coquimbo-Serena")
ciudades <- c("urb_4_18")
obs <- 1
delito <- delitos[1]

Lectura de Datos Espaciales

Para el caso del presente ejemplo práctico se utiizará la ciudad de Coquimbo-Serena.

city <- ciudades[obs]
cityp <- readRDS(paste0("data/",city,".rds"))
cityp <- spTransform(cityp, CRS(crs_utm))

El objeto espacial cityp cuyos puntos corresponde donde ha ucurrido eventos deictuales, cuyo formato SpatialPointDataFrame. A continuación se visualizan los delitos de hurto:

Matriz de vecindad espacial: Se utilizó la función propia spatial_weights() @(sec-spatial-weights)

# matriz de vecindad espacial
nb <- spatial_weights(city_df = cityp, nvec= 12)

Ventana de Trabajo para ppp (Point Pattern Plane)

Crea un objeto de clase “ppp” que representa un conjunto de datos de patrones de puntos en el plano bidimensional. Por lo anterior, se debe crear una ventana de observación w.

#make windows
  ext <- raster::extent(cityp)
  x_min <- ext[1] + 100
  x_max <- ext[2] - 100
  y_min <- ext[3] + 100
  y_max <- ext[4] - 100
  w <- as.owin(c(x_min,x_max, y_min, y_max))

Transformaciones a los datos.

Se utilizó funciónes de transformacion general mencionadas en Section 7.2 que corresponden a las siguientes:

• Transformación Cúbica
• Transformación Logarítmica
• Cáculo del Lag Espacial

cityp <- cityp %>%
    st_as_sf(crs = 4326) %>% 
    mutate(across(.cols = all_of(delitos),
                  .fns = t_cub, .names = "cub{.col}")) %>% 
    mutate(across(.cols = all_of(delitos),
                  .fns = t_log1p, .names = "log{.col}")) %>% 
    mutate(across(.cols = all_of(delitos),
                  .fns = lag_spatial, nb, .names = "lag{.col}")) %>%
    as("Spatial")

zona	id	x	y	pflot	viol	crim	abus	amen	robo	hurto	asalto	denspob	educ	tamhog	constr	ofcom	vulzon	mobil	lpflot	neduc	ldenspob	lofcom	lconstr	cubhurto	cubrobo	cubasalto	cubamen	cubabus	cubcrim	cubviol	loghurto	logrobo	logasalto	logamen	logabus	logcrim	logviol	laghurto	lagrobo	lagasalto	lagamen	lagabus	lagcrim	lagviol
4101161005	22	-71.239	-29.862	1.000	0	0	0	0	0	0	0	144.756	10.867	3.170	4328.618	0.000	0.367	0.079	0.178	0.477	0.767	0.000	0.840	0	0	0	0	0	0	0	0	0	0.000	0.000	0.000	0	0	0.083	0.083	0.083	0.417	0.250	0.083	0
4101161005	32	-71.243	-29.863	1.000	0	0	0	0	0	0	0	227.772	10.716	3.393	3830.009	0.000	0.367	0.079	0.178	0.466	0.846	0.000	0.829	0	0	0	0	0	0	0	0	0	0.000	0.000	0.000	0	0	0.000	0.083	0.083	0.917	0.583	0.083	0
4101161005	35	-71.240	-29.863	1.034	0	0	0	0	0	0	1	171.133	11.801	3.172	4328.618	0.000	0.367	0.079	0.182	0.550	0.796	0.000	0.840	0	0	1	0	0	0	0	0	0	0.693	0.000	0.000	0	0	0.000	0.083	0.000	0.500	0.250	0.083	0
4101161005	36	-71.239	-29.863	1.000	0	0	0	0	0	0	0	154.603	10.870	3.135	4328.618	0.000	0.367	0.079	0.178	0.478	0.779	0.000	0.840	0	0	0	0	0	0	0	0	0	0.000	0.000	0.000	0	0	0.083	0.083	0.083	0.417	0.250	0.083	0
4101161005	37	-71.238	-29.863	1.025	0	0	0	0	0	0	0	155.145	10.858	3.148	3230.309	0.000	0.367	0.079	0.181	0.477	0.779	0.000	0.815	0	0	0	0	0	0	0	0	0	0.000	0.000	0.000	0	0	0.083	0.083	0.083	0.000	0.250	0.083	0
4101161001	48	-71.244	-29.864	1.125	0	0	0	0	0	0	0	168.179	11.958	3.336	2384.222	14.682	0.292	0.122	0.193	0.562	0.793	0.574	0.788	0	0	0	0	0	0	0	0	0	0.000	0.000	0.000	0	0	0.000	0.000	0.250	0.500	0.417	0.083	0
4101161005	49	-71.243	-29.864	1.071	0	0	0	0	0	0	0	212.743	10.978	3.415	3795.759	0.000	0.367	0.079	0.187	0.486	0.834	0.000	0.828	0	0	0	0	0	0	0	0	0	0.000	0.000	0.000	0	0	0.000	0.083	0.250	0.917	0.583	0.083	0
4101161005	50	-71.242	-29.864	1.000	0	0	1	1	0	0	0	220.188	10.705	3.409	4026.613	0.000	0.367	0.079	0.178	0.465	0.840	0.000	0.834	0	0	0	1	1	0	0	0	0	0.000	0.693	0.693	0	0	0.000	0.083	0.000	1.000	0.417	0.083	0
4101161005	51	-71.241	-29.864	1.021	0	0	1	1	0	0	0	243.501	10.120	3.337	4328.618	0.000	0.367	0.079	0.180	0.419	0.857	0.000	0.840	0	0	0	1	1	0	0	0	0	0.000	0.693	0.693	0	0	0.000	0.083	0.083	0.667	0.250	0.000	0
4101161005	52	-71.240	-29.864	1.022	0	0	0	0	0	0	0	320.023	10.720	3.324	4328.618	0.000	0.367	0.079	0.180	0.466	0.905	0.000	0.840	0	0	0	0	0	0	0	0	0	0.000	0.000	0.000	0	0	0.000	0.083	0.167	0.583	0.250	0.000	0

4.2 Evaluación Cramer Von Mises

4.2.1 Modelos MLB

Para transformar los datos espaciales de los delitos y las predicciones de los modelos se utilizó la función propia ppp_maker() Section 7.4 que retorna el objeto ppp.

Lectura del Modelo.

model_mlb <- readRDS(paste0("output/mlb_",city,"_",delito,".rds"))

Transformaciones de las predicciones del modelo a ppp.

ppp_mbl_pred <- ppp_maker(model = model_mlb, city = cityp, 
                     delito = delito, ppp_predict = T, w = w)

plot(density(ppp_mbl_pred, adjust=.25))

Transformaciones de los delitos reales a a ppp.

ppp_mbl <- ppp_maker(model = model_mlb, city = cityp, 
                          delito = delito, ppp_predict = F, w= w)

plot(density(ppp_mbl, adjust=.25))

Evaluación cdf.test.()

cvm_mlb <- cdf.test(ppp_mbl, as.im(ppp_mbl), test="cvm")
plot(cvm_mlb)

attr_model <- attr(cvm_mlb, which = "frame")
im <- attr_model$values$Zimage 
r_mlb <- raster(im)
proj4string(r_mlb)=crs_utm
values(r_mlb) <- ifelse(values(r_mlb)==0, NA, values(r_mlb))
plot(r_mlb)

pred_mlb <- predict_df(model = model_mlb, city = cityp, longlat = F)
pred_mlb <- st_as_sf(x = pred_mlb,                         
                       coords = c("x", "y"),
                       crs = 32719) %>% 
  filter(!is.na(delitos))

mapview(r_mlb, na.color= NA)+
  mapview(pred_mlb, zcol = "ypred",hide = TRUE,  cex = 2)+
  mapview(pred_mlb, zcol = "delitos", hide = TRUE,  cex = 2)

Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
+b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
+wktext +no_defs +type=crs

Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

4.2.2 Modelos SAR

city <- ciudades[obs] 
model_sar <- readRDS(paste0("output/sar_",city,"_",delito,".rds"))

ppp_sar_pred <- ppp_maker(model = model_sar, city = cityp, 
                            delito = delito, ppp_predict = T,
                            w = w, nb = nb)

This method assumes the response is known - see manual page

plot(density(ppp_sar_pred, adjust=.1))

ppp_sar <- ppp_maker(model = model_sar, city = cityp, 
                       delito = delito, ppp_predict = F, 
                       w = w, nb = nb)
  
plot(density(ppp_sar, adjust=.25))

# st <- syrjala.test(ppp_mbl, ppp_mbl_pred, nsim = 999)
cvm_sar <- cdf.test(ppp_sar, as.im(ppp_sar_pred), test="cvm")
plot(cvm_sar)

attr_model <- attr(cvm_sar, which = "frame")
im <- attr_model$values$Zimage 
r_sar <- raster(im)
proj4string(r_sar)=crs_utm
values(r_sar) <- ifelse(values(r_sar)==0, NA, values(r_sar))
# plot(r_sar)

pred_sar <- predict_df(model = model_sar, city = cityp, longlat = F)

This method assumes the response is known - see manual page

pred_sar <- st_as_sf(x = pred_sar,                         
                       coords = c("x", "y"),
                       crs = 32719) %>% 
  filter(!is.na(delitos))

mapview(r_sar, na.color= NA)+
  mapview(pred_sar, zcol = "ypred",hide = TRUE,  cex = 2)+
  mapview(pred_sar, zcol = "delitos", hide = TRUE,  cex = 2)